﻿Imports WINX32Lib 'Remember to add r3eference first (Add Reference option). For MS Visual Studio .NET 2003: before adding a reference to the project, corresponding assembly should be added to the Add Reference dialog box (http://msdn.microsoft.com/en-us/library/wkze6zky%28v=vs.71%29.aspx).
Imports System.IO
Imports System.Math


Module Grain_nalysis

    Sub Analyze_Frame_Grains(ByRef Frame_To_Analyze As Array, ByRef Frame_DimX As Integer, ByRef Frame_DimY As Integer, _
                     ByRef Frame_Index As Integer, _
                 ByVal low As Short, ByVal high As Short, ByVal Min_Intensity As Short, ByVal low1 As Short, ByVal high1 As Short, _
                 Optional ByRef Spectrum(,) As Single = Nothing, _
                 Optional ByRef Good_Events_Sum()(,) As Single = Nothing, _
                 Optional ByRef Total_Number_Of_Good_Events(,) As Integer = Nothing, _
                 Optional ByRef Intensity_Histograme()() As Integer = Nothing, _
                 Optional ByRef Multiple_Hits()(,) As Integer = Nothing, _
                 Optional ByRef Total_Number_Of_Multiple_Hits(,) As Integer = Nothing, Optional ByRef toc2 As String = "")

        Dim tic As String
        tic = (Now().Millisecond) / 1000 + Now().Second + (Now().Minute) * 60

        If Not Frame_To_Analyze.Rank = 2 Then
            MsgBox("The frame to analize hase a dimention of " & Frame_To_Analyze.Rank)
        End If

        Dim Max_Intensity As Integer
        Dim Max_Multiplicity As Integer = low1 \ (high1 - low1)
        Dim Neighbours_Sum As Integer
        Dim Frame_As_Grain_Numbers(Frame_DimX, Frame_DimX) As Integer

        'We need to clear out the edges (alternatively we could put our array into a bigger one with zeros around)
        For i = 0 To Frame_DimX
            Frame_To_Analyze(i, 0) = 0
            Frame_To_Analyze(i, Frame_DimY) = 0
        Next
        For j = 1 To Frame_DimY - 1
            Frame_To_Analyze(0, j) = 0
            Frame_To_Analyze(Frame_DimX, 0) = 0
        Next

        ' Number of elements in Frame (also used as maximum (never reached) grain)
        Dim NumEl = 1000 'Frame_DimX * Frame_DimY - 1

        'Here we check wich pixels should be analize -
        ' we choose only the ones that have the intensity abowe the treshold 
        ' and with  the sum of naighbourf abowe low1er limit of the energy window 
        For i = 1 To Frame_DimX - 1
            For j = 1 To Frame_DimY - 1

                'Intensity_Histograme(0)(Frame_To_Analyze(i, j)) += 1
                If Frame_To_Analyze(i, j) > Min_Intensity Then
                    Neighbours_Sum = Frame_To_Analyze(i - 1, j - 1) + _
                     Frame_To_Analyze(i - 1, j) + _
                     Frame_To_Analyze(i - 1, j + 1) + _
                     Frame_To_Analyze(i, j - 1) + _
                     Frame_To_Analyze(i, j) + _
                     Frame_To_Analyze(i, j + 1) + _
                     Frame_To_Analyze(i + 1, j - 1) + _
                     Frame_To_Analyze(i + 1, j) + _
                     Frame_To_Analyze(i + 1, j + 1)
                    Intensity_Histograme(2)(Neighbours_Sum) += 1
                    If Neighbours_Sum > low1 Then
                        Dim ii, jj As Integer ' this will be the indexing for Frame_As_Grain_Numbers - which is bigger
                        ii = i
                        jj = j
                        Frame_As_Grain_Numbers(ii - 1, jj - 1) = NumEl
                        Frame_As_Grain_Numbers(ii - 1, jj) = NumEl
                        Frame_As_Grain_Numbers(ii - 1, jj + 1) = NumEl
                        Frame_As_Grain_Numbers(ii, jj - 1) = NumEl
                        Frame_As_Grain_Numbers(ii, jj) = NumEl
                        Frame_As_Grain_Numbers(ii, jj + 1) = NumEl
                        Frame_As_Grain_Numbers(ii + 1, jj - 1) = NumEl
                        Frame_As_Grain_Numbers(ii + 1, jj) = NumEl
                        Frame_As_Grain_Numbers(ii + 1, j + 1) = NumEl
                    End If
                End If
            Next
        Next

        Dim actual_grain_number As Integer = 1
        Dim grain_intensity() As Integer
        Dim number_of_pixels_in_grain() As Integer
        Dim Grain_Mutiplicity(,) As Integer
        ReDim grain_intensity(NumEl)
        ReDim number_of_pixels_in_grain(NumEl)
        ReDim Grain_Mutiplicity(NumEl, 2)
        For i = 0 To Frame_DimX
            For j = 0 To Frame_DimY
                If Frame_As_Grain_Numbers(i, j) > actual_grain_number Then
                    ' Sum of intensities:
                    Dim Sum = 0
                    ' Start with the original point in the stack.
                    Dim pts As New Stack(Of Integer())(1000)
                    Dim pixel() As Integer = {i, j}
                    pts.Push(pixel)
                    Frame_As_Grain_Numbers(i, j) = actual_grain_number
                    ' While the stack is not empty, process a pixel.
                    Do While pts.Count > 0
                        Dim pt() As Integer = pts.Pop()
                        Sum += Frame_To_Analyze(pt(0), pt(1))
                        number_of_pixels_in_grain(actual_grain_number) += 1
                        If pt(0) > 0 Then SafeCheckPixel(Frame_As_Grain_Numbers, pts, pt(0) - _
                            1, pt(1), NumEl, actual_grain_number)
                        If pt(1) > 0 Then SafeCheckPixel(Frame_As_Grain_Numbers, pts, pt(0), _
                            pt(1) - 1, NumEl, actual_grain_number)
                        If pt(0) < Frame_As_Grain_Numbers.Length - 1 Then SafeCheckPixel(Frame_As_Grain_Numbers, _
                            pts, pt(0) + 1, pt(1), NumEl, actual_grain_number)
                        If pt(1) < Frame_As_Grain_Numbers.Length - 1 Then SafeCheckPixel(Frame_As_Grain_Numbers, _
                            pts, pt(0), pt(1) + 1, NumEl, actual_grain_number)
                    Loop

                    grain_intensity(actual_grain_number) = Sum
                    '#################### end of algorithm #########
                    Max_Intensity = Max(Max_Intensity, Sum)
                    Intensity_Histograme(1)(Sum) += 1

                    Dim Lower_Lewel_Multiplicity As Integer = Sum \ low1
                    Dim Upper_Lewel_Multiplicity As Integer = Ceiling(Sum / high1)
                    'here we check either we are in an energy window
                    If Lower_Lewel_Multiplicity >= Upper_Lewel_Multiplicity Then
                        Total_Number_Of_Good_Events(Frame_Index, 0) += Lower_Lewel_Multiplicity
                        Grain_Mutiplicity(actual_grain_number, 0) = Lower_Lewel_Multiplicity
                        If Lower_Lewel_Multiplicity = Upper_Lewel_Multiplicity Then
                            Total_Number_Of_Multiple_Hits(Frame_Index, Lower_Lewel_Multiplicity) += 1
                        Else
                            Total_Number_Of_Good_Events(Frame_Index, 1) += Lower_Lewel_Multiplicity - Upper_Lewel_Multiplicity
                            Grain_Mutiplicity(actual_grain_number, 1) = Lower_Lewel_Multiplicity - Upper_Lewel_Multiplicity
                        End If
                    End If
                    actual_grain_number += 1
                End If
            Next
        Next

        'Corrected frame And Spectrum
        For i = 0 To Frame_DimX
            For j = 0 To Frame_DimY
                If Frame_To_Analyze(i, j) > 0 Then
                    Intensity_Histograme(0)(Frame_To_Analyze(i, j)) += 1
                End If
                If Frame_As_Grain_Numbers(i, j) > 0 Then
                    Good_Events_Sum(0)(i, j) += _
                        Grain_Mutiplicity(Frame_As_Grain_Numbers(i, j), 0) / number_of_pixels_in_grain(Frame_As_Grain_Numbers(i, j))
                    Good_Events_Sum(1)(i, j) += _
                        Grain_Mutiplicity(Frame_As_Grain_Numbers(i, j), 1) / number_of_pixels_in_grain(Frame_As_Grain_Numbers(i, j))
                    Spectrum(i, 0) += Grain_Mutiplicity(Frame_As_Grain_Numbers(i, j), 0) / number_of_pixels_in_grain(Frame_As_Grain_Numbers(i, j))
                    Spectrum(i, 1) += Grain_Mutiplicity(Frame_As_Grain_Numbers(i, j), 1) / number_of_pixels_in_grain(Frame_As_Grain_Numbers(i, j))
                ElseIf Frame_To_Analyze(i, j) > 0 Then
                    Intensity_Histograme(1)(Frame_To_Analyze(i, j)) += 1
                End If
            Next
        Next

        toc2 = (Now().Millisecond) / 1000 + Now().Second + (Now().Minute) * 60 - tic
    End Sub


    ' See if this point should be added to the stack.
    Private Sub SafeCheckPixel(ByVal bm(,) As Integer, ByVal pts As  _
        Stack(Of Integer()), ByVal x As Integer, ByVal y As Integer, ByVal _
        old_value As Integer, ByVal new_value As Integer)
        Dim clr As Integer = bm(x, y)
        If clr.Equals(old_value) Then
            Dim point() As Integer = {x, y}
            pts.Push(point)
            bm(x, y) = new_value
        End If
    End Sub


End Module

